Why a talk on programming language ?¶
- Better Software, Better Research (Goble, C. (2014). IEEE Internet Computing)
Of 2,000 scientists surveyed online, 91% said using scientific software is important for their own research, 84% said developing scientific software is important for their own research [...] and 38% spend at least one fifth of their time developing software.
“You can download our code from the URL supplied. Good luck downloading the only postdoc who can get it to run, though” (Ian Holmes)
- Big Data -> Big Software
- Loi pour une république numérique (2017) $\rightarrow$ new codes must (more or less) be open source (cf https://scinfolex.com/2017/12/08/les-logiciels-produits-par-les-administrations-sont-passes-en-open-source-par-defaut-et-voici-pourquoi/)
Outline¶
- State of the art : the two language problem
- What is Julia
- Why Julia is fast
- Julia Ecosystem & interop
- (Parallelism)
State of the art : Matlab¶
Pros :
- Polished product with support
- Simulink
- High level syntax
Cons :
- Closed : not everyone has access to it, impossible to put a demonstrator online.
- Slow loops (better since 2015) : not everything is pretty once vectorized
- Not fast per se (Fortran bindings), memory management is difficult
- Not a generalist language (cumbersome to put a demo on the web e.g.)
- 12 000€ each year for the lab
State of the art : Python¶
Pros :
- Free and open-source, generalist, widly used outside scientific community
- Lot of scientific communities are embracing it
- Lots of efforts to make it fast (numba, ...)
Cons :
- Scientific computing is not native :
- all fallback to C/Fortran black-boxes -> limit flexibility
- Object Oriented paradigm can be cumbersome for scientific code
Science and the two languages problem¶
Scientists need to easily explore new ideas :
- Need for mathematical abstractions
Need for customizations (avoid black boxes, try variations)
But also need for performance (intrinsic goal or to allow re-use)
What is done now (when you need to go further than using existing packages, for e.g. data analysis)
1) Prototyping in R/Python/Matlab
2) Rewriting whole or parts (if money/man power) to a performant low-level language as C/Fortran
Yet another language¶
Here comes Julia :¶
- innovative open-source generalist programming language with scientific computing at its core
- easy as Matlab, fast as Fortran, flexible as Python, deep as LISP
- leverages (for now...) major C/Fortran science packages, e.g. LAPACK, MPI, FFTW...
- 5th language to run at 1 petaflops (Celeste project), after assembly, Fortran, C, C++
- State of the art packages in Optimization (JUMP), Differential Equations (DifEq), ... Well positioned for ML (Flux, Knet, autodiff...)
- solves the "two-language problem"
Background¶
- origins at MIT, UCSB applied math & computer science, 2010
- founders Viral Shah, Jeff Bezanson, Stefan Karpinsky (now at Julia Computing LLC), Alan Edelman (MIT)
- ~10-person core language team, ~600 contributors, ~1500 registered packages
- support from Intel, Microsoft, Wall St., Moore Foundation, NumFocus
- julia-0.1 released in 2012, julia-1.0 in August 2018
Its goal : Solve the languages problem by having a dynamic, high level language with good mathematical expressivity able to be as fast as C/Fortran.
Disclaimer¶
- I am not a Julia expert : I've been following the language for ~ one year and started porting Python codes a few months ago.
- Julia (like Python...) is not a drop-in replacement for Matlab, it is an alternative.
- Don't fix that ain't broken :
- if Matlab/R/Python/C is good for you, no need to change
- but if you are sometimes faced with some limitations of these languages, give Julia a try
In [1]:
using LinearAlgebra
A = randn(4,4) # 4 x 4 matrix of normally distributed random numbers
Out[1]:
In [2]:
A[:,1] # familiar colon syntax --extract first column of A
Out[2]:
Linear Algebra¶
In [3]:
b = randn(4)
x = A\b # solve Ax=b for x using backslash operator
norm(A*x-b)
Out[3]:
In [4]:
U, σ, V = svd(A); # singular value decomposition --note unicode variable name σ
In [5]:
Σ = diagm(0 => σ)
norm(A - U*Σ*V') # compute error of SVD factorization A = U Σ V'
Out[5]:
Fast as Fortran¶
In [6]:
function loop()
a=0
for i in 1:1_000_000
a+=1;
end
end
@time loop # compile time
@time loop # now fast
Out[6]:
Compare with Matlab (R2018)¶
function [] = loop()
a=0;
for i=1:1000000
a=i+1;
end
end
f = @() loop();
timeit(f) => 0.0018s
2000x slower !
Benchmarks on common code pattterns¶
Takeaway: Julia timings are clustered around C timing, Matlab/Python/R orders of magnitude slower.
KS-CNAB2 benchmark: CPU time versus lines of code¶
Takeaway: The Julia PDE code is almost as compact as Matlab/Python, almost as fast as C/Fortran.
Julia: easy, dynamic, and fast. How?¶
- Just-in-time compilation (JIT)
- user-level code is compiled to machine code on-the-fly
- Meticulous type system
- designed to maximize impact of JIT
- type inference: compiler determines types of variables
- Multiple dispatch
- functions are compiled for each set of argument types
- function dispatch determined at compile time when possible, run time when not
Just-in-time compilation¶
Functions are compiled to machine code when first run. Subsequent runs are as fast as compiled C, C++, Fortran.
Caveats¶
- It's a marathon not a sprint : Julia can be slow at first
- It's an open-source language with lots of good willing contributors, it's not a product (see JuliaPro...)
- If you directly translate from Matlab to Julia you will not always see an improvement, you have exploit Julia strengths (loops to avoid allocations, type stability, ...)
Logistic map example¶
$$f(x) = 4\, x\, (1-x)$$In [7]:
f(x) = 4x*(1-x) # define logistic map
@time f(0.3); # first run is slow
@time f(0.4); # second run is a thousand times faster
Observe the compilation of $f(x)$ by stages¶
user Julia -> generic Julia expression -> typed Julia expression -> intermediate compiler language -> machine code
In [8]:
include("macros.jl");
In [9]:
@code_lowered f(0.3) # show f(x) as generic Julia expression
Out[9]:
In [10]:
@code_typed f(0.3) # show f(x) as Julia expr with inferred types, based on the arg types
Out[10]:
Code transformation¶
The code is then converted to intermediate compiler language (LLVM), and finally compiled by LLVM to machine code
In [11]:
@code_llvm_nocomment f(0.3) # show f(x) in intermediate compiler language (LLVM)
In [12]:
@code_native_nocomment f(0.3) # show f(x) in IA-64 assembly language
Type inference and dispatch¶
We can see that if we call f on an Integer we get a code specialised for integer (i64)
In [11]:
@code_llvm_nocomment f(0.3) # show f(x) in intermediate compiler language (LLVM)
In [13]:
@code_llvm_nocomment f(3) # show f(x) in intermediate compiler language (LLVM)
Type stability¶
No need to declare type but need to be type stable...
In [14]:
function g1(a,n)
x = 0
for i in 1:n
x += a
end
end
@code_warntype g1(2.0,10) # see the Union{Float64, Int64}
# @code_llvm_nocomment g1(2.0,10)
In [15]:
function g2(a,n)
x = zero(a) # get the zero of the type of a
for i in 1:n
x += a
end
end
@code_warntype g2(2.0,10) # the type stable version is shorter thus faster
#@code_llvm_nocomment g2(2.0,10)
Millionth iterate of logistic map: comparison of Julia, C, Matlab¶
Given $f(x) = 4\,x \, (1-x)$,
let $~~f^N(x) = f(\,f(\,...(\,f(\,f(x)))))$ be the millionth iterate of $f$, with $N=10^6$.
$f^N$ in Julia: run-time code generation & compilation!¶
In [16]:
# define function that, given function g(x),
# returns iterated function gᴺ(x) = g(g(...(g(g(x)))))
function iterator(g, N)
# construct gᴺ, the Nth iterate of g
function gᴺ(x)
for i ∈ 1:N
x = g(x)
end
return x
end
return gᴺ
end
fᴺ = iterator(f, 10^6); # generate millionth iterate funcion fᴺ(x) for f(x) = 4x(1-x)
fᴺ(0.67) # run fᴺ(x) once to compile it
Out[16]:
$f^N$ in C++, compiled
In [17]:
; cat codes/fN.cpp
In [18]:
print("t=");
@time x = fᴺ(0.67); # run the Julia fᴺ on x = 0.67
@show x;
C++¶
In [19]:
; g++ -O3 -o fN codes/fN.cpp -lm
In [20]:
; ./fN 0.67
Julia and C++ get the same result x in roughly same execution time t. Sometimes Julia is faster, sometimes C++.
Matlab¶
>> tic(); x=fN(0.67); t=toc();
>> t, x
t = 0.048889000000000
x = 0.101168853345475Matlab gets the same result $x$, but ten to twenty times slower than Julia or C++.
A short break¶
An ecosystem of packages¶
- Most of packages are available on Github (ease collaboration)
- Main fields are grouped in Github Organizations (see https://julialang.org/community/)
- Julia comes with a powerful package/envrionment manager
In [21]:
using Pkg #load the Pkg stdlib
Pkg.activate(".") #activate the local environment
Out[21]:
Singular values of Hilbert matrix and arbitrary precision¶
The Hilbert matrix $A_{ij} = (i+j-1)^{-1}$ is notoriously ill-conditioned.
In [22]:
using LinearAlgebra #load stdlib
#install package : no need online, needed on your computer
#Pkg.add(PackageSpec(name="GenericSVD",rev="master"));
using GenericSVD #load package
In [23]:
m = 8
A = [1//(i+j-1) for i=1:m, j=1:m] # 8 x 8 Hilbert matrix of Rationals
Out[23]:
In [24]:
# 16 x 16 Hilbert is too ill-conditioned for 64-bit arithmetic
m = 16 ; A = [1/(i+j-1) for i=1:m, j=1:m] # 16 x 16 Hilbert matrix of Float64s
@show cond(A)
@show eps(Float64);
In [25]:
σ = svdvals(A)
σ[end-5:end] # last 5 values
Out[25]:
In [26]:
# Make 32 x 32 Hilbert matrix of 256-bit BigFloats and show a few elements
setprecision(256)
m = 32
A = [BigFloat(1//(i+j-1)) for i=1:m, j=1:m];
@show A[1,1]
@show A[2,1]
@show A[3,1]
;
In [27]:
# Compute singular values of 32 x 32 Hilbert matrix in 256-bit arithmetic
σ = svdvals!(A);
σ[end-5:end] # last 5 values
Out[27]:
User-defined type: Galois field GF(p)¶
Finite scalar field GF(p) is $\mathbb{Z}/p$, the integers modulo p, where p is prime. Example by Andreas Noack, Julia Computing.
In [28]:
# Type definition: Galois fields GF(p), where p is prime modulus, T is integer type
struct GF{p,T} <: Number where {p,T<:Integer}
rep::T # representative integer which holds the value of a GF(p) variable
function GF{p,T}(x::Integer) where {p,T<:Integer}
return new(mod(x, p))
end
end
GF{p}(x::T) where {p,T<:Integer} = GF{p,T}(x)
In [48]:
# Define some basic methods for GF(p) that all Julia Numbers must have
import Base: convert, inv, one, promote_rule, show, zero
function call(::Type{GF{p}}, x::Integer) where {p}
if !isprime(p)
throw(ArgumentError("p must be prime"))
end
return GF{p,typeof(x)}(mod(x, p))
end
convert(::Type{GF{p,T}}, x::Integer) where {p,T} = GF{p}(x)
convert(::Type{GF{p}}, x::Integer) where {p} = GF{p}(x)
convert(::Type{GF{p,T}}, x::GF{p}) where {p,T} = GF{p,T}(x.rep)
promote_rule(::Type{GF{p,T1}}, ::Type{T2}) where {p,T1,T2<:Integer} = GF{p,promote_type(T1,T2)}
show(io::IO, x::GF) = show(io, x.rep);
In [29]:
# Define arithmetic operations on GF(p)
import Base: +, -, *, /,conj
for op in (:+, :-, :*) # loop over ops, defining each in terms of integer ops mod p
@eval begin
($op)(x::GF{p,T}, y::GF{p,T}) where {p,T} = GF{p,T}($(op)(x.rep, y.rep))
end
end
# Define inverse and ÷. Requires more care than +, -, * because dividing by zero throws error
function inv(x::GF{p,T}) where {p,T}
if x == zero(x)
throw(DivideError())
end
r, u, v = gcdx(x.rep, p)
GF{p,T}(u)
end
function conj(x::GF{p,T}) where {p,T}
u = conj(x.rep)
GF{p,T}(u)
end
(/)(x::GF{p}, y::GF{p}) where {p}= x*inv(y);
In [30]:
# Create some GF(7) variables and do arithmetic
x = GF{7}(9) # x = 9 mod 7 = 2
y = GF{7}(12) # y = 12 mod 7 = 5
@show x
@show y
@show x + y # 2 + 5 mod 7 = 0
@show x - y # 2 - 5 mod 7 = 4
@show x * y # 2 * 5 mod 7 = 3
@show x / y # 2 ÷ 5 mod 7 = 6, because 2 = 6*5 mod 7
;
Linear algebra over Galois field¶
In [31]:
import Random: seed!
seed!(1234)
A = [GF{7}(rand(0:6)) for i = 1:4, j = 1:4]; # matrix of random GF(7) elems
In [32]:
b = [GF{7}(rand(0:6)) for i = 1:4]; # random vector b for Ax=b problem
In [33]:
x = A\b; # solve Ax=b over GF(p)!
In [34]:
A*x - b # check that x satisfies Ax=b
Out[34]:
Whoa! Built-in backslash operator on matrix of GF(p) worked!
Julia generated and compiled a GF(p)-specific version of its generic LU decomp function. All it needs is definitions of $+, -, \times, \div$.
Solving Differential Equations¶
Solve the Lorenz equations:
$$ \begin{align} \frac{dx}{dt} &= σ(y-x) \\ \frac{dy}{dt} &= x(ρ-z) - y \\ \frac{dz}{dt} &= xy - βz \\ \end{align} $$In [35]:
#Pkg.add("DifferentialEquations") # add the Differential equation suite
using DifferentialEquations # first time is very slow (precompilation)
using Plots